Classification of Variants across paired longitudinal samples in the PBTA Cohort

Set up

suppressPackageStartupMessages({
  library(tidyverse)
})

Directories and File Inputs/Outputs

# Detect the ".git" folder -- this will be in the project root directory
# Use this as the root directory to ensure proper sourcing of functions 
# no matter where this is called from
root_dir <- rprojroot::find_root(rprojroot::has_dir(".git"))
analysis_dir <- file.path(root_dir, "analyses", "tmb-vaf-longitudinal")
results_dir <- file.path(analysis_dir, "results")
input_dir <- file.path(analysis_dir, "input")
files_dir <- file.path(root_dir, "analyses", "sample-distribution-analysis", "results")

# Input files
genomic_paired_file <- file.path(files_dir, "genomic_assays_matched_time_points_extended.tsv") 
tmb_vaf_file <- file.path(results_dir, "tmb_vaf_genomic.tsv")
palette_file <- file.path(root_dir, "figures", "palettes", "oncoprint_color_palette.tsv")

# File path to plot directory
plots_dir <-
  file.path(analysis_dir, "plots")
if (!dir.exists(plots_dir)) {
  dir.create(plots_dir)
}

source(paste0(root_dir, "/figures/scripts/theme.R"))

Read in data and process

tmb_vaf_df <- readr::read_tsv(tmb_vaf_file, guess_max = 100000, show_col_types = FALSE) %>% 
  filter(!tmb >= 10) %>% 
  select(Kids_First_Biospecimen_ID, Variant_Classification, gene_protein, mutation_count,   region_size, tmb, VAF)

genomic_paired_df <- readr::read_tsv(genomic_paired_file, guess_max = 100000, show_col_types = FALSE) %>% 
  left_join(tmb_vaf_df, by = c(\Kids_First_Biospecimen_ID\)) %>%
  filter(!is.na(tmb)) 


# Attention as some bs specimen might not have TMB!
# If that happens, we will end up with samples lacking timepoints.

# Which patient samples don't have TMB?
# genomic_paired_df %>% 
#  filter(is.na(tmb)) %>% 
#  unique() %>% 
#  regulartable() %>%
#  fontsize(size = 12, part = \all\)

descriptors_df <- genomic_paired_df %>%
  group_by(Kids_First_Participant_ID) %>%
  summarize(descriptors = paste(sort(tumor_descriptor), collapse = \

# Read color palette
palette_df <- readr::read_tsv(palette_file, guess_max = 100000, show_col_types = FALSE) 

# Define and order palette
palette <- palette_df$hex_codes
names(palette) <- palette_df$Variant_Classification
# Define label for plots
Alteration_type <- df_plot$Variant_Classification

# Define ylim
ylim <- max(df_plot$log10_tmb)

What type of alterations we observe per tumor descriptor?

# Create bxp
print(ggpubr::ggboxplot(df_plot, 
                        x = \tumor_descriptor\, 
                        y = \log10_tmb\, 
                        color = \Variant_Classification\,
                        palette = palette) +
        theme_Publication() + 
        scale_y_continuous(limits = c(0, ylim)) +
        xlab(\Timepoint\) +
        theme(axis.text.x = element_text(angle = 90)))

# Save the plot
ggsave(filename = \Alteration_type_timepoints.pdf\, 
       path = plots_dir, 
       width = 15, 
       height = 8, 
       device = \pdf\, 
       useDingbats = FALSE)

What type of alterations we observe per tumor descriptor in each cancer group?

# Create bxp
print(ggpubr::ggboxplot(df_plot, 
                        x = \tumor_descriptor\, 
                        y = \log10_tmb\, 
                        color = \Variant_Classification\,
                        palette = palette) +
        facet_wrap(~cg_id) +
        theme_Publication() + 
        xlab(\Timepoint\) +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90)))

# Save the plot
ggsave(filename = \Alteration_type_cancer_group.pdf\, 
       path = plots_dir, 
       width = 25, 
       height = 18, 
       device = \pdf\, 
       useDingbats = FALSE)

What type of alterations we observe per tumor descriptor in each cancer group defined by cgGFAC?

df_plot_cgGFAC <- df_plot %>% 
  arrange(tumor_descriptor_cgGFAC_n)
  #mutate(tumor_descriptor_cgGFAC_n = factor(tumor_descriptor_cgGFAC_n)) 

#df_plot_cgGFAC$tumor_descriptor_cgGFAC_n %>% levels()

# Create bxp
print(ggpubr::ggboxplot(df_plot_cgGFAC, 
                        x = \tumor_descriptor_cgGFAC_n\, 
                        y = \log10_tmb\, 
                        color = \Variant_Classification\,
                        palette = palette) +
        theme_Publication() + 
        xlab(\Timepoint\) +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90)))

# Save the plot
ggsave(filename = \Alteration_type_cgGFAC.pdf\, 
       path = plots_dir, 
       width = 14, 
       height = 8, 
       device = \pdf\, 
       useDingbats = FALSE)

cgGFAC_id <- as.character(unique(df_plot_cgGFAC$cgGFAC))
cgGFAC_id

[1] \ATRT\  \DMG\   \HGG\   \LGG\   \Other\
ATRT

DMG

HGG

LGG

Other

# Loop through variable
for (i in seq_along(cgGFAC_id)){
  print(i)
  df_sub <- df_plot_cgGFAC %>%
      filter(cgGFAC == cgGFAC_id[i])

  
   # Create bxp
  print(ggpubr::ggboxplot(df_sub, 
                        x = \tumor_descriptor_cgGFAC_n\, 
                        y = \log10_tmb\, 
                        color = \Variant_Classification\,
                        palette = palette) +
        theme_Publication() + 
        labs(title = paste(cgGFAC_id[i])) +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)))
}

[1] 1

[1] 2

[1] 3

[1] 4

[1] 5

What type of alterations we observe per tumor descriptor in each cancer group (add _n))?

cg <- as.character(unique(df_plot$cg_id))
cg

 [1] \Embryonal tumor with multilayer rosettes\
 [2] \Low-grade glioma\                        
 [3] \Ependymoma\                              
 [4] \Medulloblastoma\                         
 [5] \Atypical Teratoid Rhabdoid Tumor\        
 [6] \Diffuse midline glioma\                  
 [7] \High-grade glioma\                       
 [8] \Ganglioglioma\                           
 [9] \Meningioma\                              
[10] \Pilocytic astrocytoma\                   
[11] \CNS Embryonal tumor\                     
[12] \Neuroblastoma\                           
[13] \Schwannoma\                              
[14] \Chordoma\                                
[15] \Malignant peripheral nerve sheath tumor\ 
[16] \Choroid plexus carcinoma\                
[17] \Adamantinomatous Craniopharyngioma\      
[18] \Dysembryoplastic neuroepithelial tumor\  
[19] \Ewing sarcoma\                           
[20] \Rosai-Dorfman disease\                   
[21] \Neurofibroma/Plexiform\                  
[22] \Glial-neuronal tumor\                    
[23] \Hemangioblastoma\                        
[24] \Craniopharyngioma\                       
Embryonal tumor with multilayer rosettes

Low-grade glioma

Ependymoma

Medulloblastoma

Atypical Teratoid Rhabdoid Tumor

Diffuse midline glioma

High-grade glioma

Ganglioglioma

Meningioma

Pilocytic astrocytoma

CNS Embryonal tumor

Neuroblastoma

Schwannoma

Chordoma

Malignant peripheral nerve sheath tumor

Choroid plexus carcinoma

Adamantinomatous Craniopharyngioma

Dysembryoplastic neuroepithelial tumor

Ewing sarcoma

Rosai-Dorfman disease

Neurofibroma/Plexiform

Glial-neuronal tumor

Hemangioblastoma

Craniopharyngioma

# Loop through variable
for (i in seq_along(cg)){
  print(i)
  df_sub <- df_plot %>%
      filter(cg_id == cg[i])
  
  # Create bxp
  print(ggpubr::ggboxplot(df_sub, 
                        x = \tumor_descriptor_cg_n\, 
                        y = \log10_tmb\, 
                        color = \Variant_Classification\,
                        palette = palette) +
        theme_Publication() + 
        xlab(\Timepoint\) +
        labs(title = paste(cg[i])) +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)))
}

[1] 1

[1] 2

[1] 3

[1] 4

[1] 5

[1] 6

[1] 7

[1] 8

[1] 9

[1] 10

[1] 11

[1] 12

[1] 13

[1] 14

[1] 15

[1] 16

[1] 17

[1] 18

[1] 19

[1] 20

[1] 21

[1] 22

[1] 23

[1] 24

What type of alterations we observe per tumor descriptor in each cancer group and timepoint model?

tm <- as.character(unique(df_plot$timepoints_models))
tm

 [1] \Dx-Rec\         \Dx-Pro\         \Dx-Dec\         \Pro-Dec\       
 [5] \Rec-SM\         \Pro-Rec\        \Dx-SM\          \Pro-Rec-Dec\   
 [9] \Dx-Pro-Rec\     \Rec-Dec\        \Dx-Pro-Rec-Dec\ \Dx-Rec-Dec\    
[13] \Dx-Pro-Dec\    
Dx-Rec

Dx-Pro

Dx-Dec

Pro-Dec

Rec-SM

Pro-Rec

Dx-SM

Pro-Rec-Dec

Dx-Pro-Rec

Rec-Dec

Dx-Pro-Rec-Dec

Dx-Rec-Dec

Dx-Pro-Dec

# Loop through variable
for (i in seq_along(tm)){
  print(i)
  df_sub <- df_plot %>%
      filter(timepoints_models == tm[i])
  
  # Create bxp
  print(ggpubr::ggboxplot(df_sub, 
                        x = \tumor_descriptor\, 
                        y = \tmb\, 
                        color = \Variant_Classification\,
                        palette = palette) +
        facet_wrap(~cancer_group) +
        theme_Publication() +
        ylab(\TMB\) +
        xlab(\Timepoint\) +
        scale_y_continuous(limits = c(0, ylim)) +
        theme(axis.text.x = element_text(angle = 90)))
}

[1] 1

[1] 2

[1] 3

[1] 4

[1] 5

[1] 6

[1] 7

[1] 8

[1] 9

[1] 10

[1] 11

[1] 12

[1] 13

sessionInfo()

R version 4.2.3 (2023-03-15)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.5.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggthemes_4.2.4  lubridate_1.9.2 forcats_1.0.0   stringr_1.5.0  
 [5] dplyr_1.1.2     purrr_1.0.1     readr_2.1.4     tidyr_1.3.0    
 [9] tibble_3.2.1    ggplot2_3.4.2   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.0  xfun_0.39         bslib_0.5.0       carData_3.0-5    
 [5] colorspace_2.1-0  vctrs_0.6.3       generics_0.1.3    htmltools_0.5.5  
 [9] yaml_2.3.7        utf8_1.2.3        rlang_1.1.1       jquerylib_0.1.4  
[13] pillar_1.9.0      ggpubr_0.6.0      glue_1.6.2        withr_2.5.0      
[17] bit64_4.0.5       lifecycle_1.0.3   munsell_0.5.0     ggsignif_0.6.4   
[21] gtable_0.3.3      ragg_1.2.5        evaluate_0.21     labeling_0.4.2   
[25] knitr_1.43        tzdb_0.4.0        fastmap_1.1.1     parallel_4.2.3   
[29] fansi_1.0.4       highr_0.10        broom_1.0.5       scales_1.2.1     
[33] backports_1.4.1   cachem_1.0.8      vroom_1.6.3       jsonlite_1.8.7   
[37] abind_1.4-5       systemfonts_1.0.4 farver_2.1.1      bit_4.0.5        
[41] textshaping_0.3.6 hms_1.1.3         digest_0.6.33     stringi_1.7.12   
[45] rstatix_0.7.2     rprojroot_2.0.3   cli_3.6.1         tools_4.2.3      
[49] magrittr_2.0.3    sass_0.4.7        crayon_1.5.2      car_3.1-2        
[53] pkgconfig_2.0.3   timechange_0.2.0  rmarkdown_2.23    R6_2.5.1         
[57] compiler_4.2.3